clear
clear all
cap log close

log using "log-files/2_Figure_1.smcl", replace smcl

timer on 2

*----------*
* Figure 1 *
*----------*
	* Figure 1a
	*----------
		* Master data
		use "data/outputs/computo_ulttranstrepdate_missingsatmunmedian.dta", clear

		* Trimming
		cumul trep_date_ulttrans, gen(cum_date)
		drop if cum_date <= 0.02 | cum_date >= 0.98
		drop cum_date
			
		* Margin
		gen mmargin_nbnn_so = mshare_nbnn_so - cshare_nbnn_so

		* Without this, Stata can't calculate the right ROT bandwidth
		su trep_date_ulttrans, d
		local global_mean = r(mean)
		local global_sd = r(sd)
		gen trep_date_ulttrans_norm = (trep_date_ulttrans - `global_mean') / `global_sd'

		* Restrict to regression sample
		quietly reg mmargin_nbnn_so $controls_educ trep_date_ulttrans_norm
		keep if e(sample)
		
		* Binned data using binsreg and all booths
		binsreg mmargin_nbnn_so trep_date_ulttrans_norm, ///
			savedata(temp) replace

		* Tack binned data output on to the end of the main data set	
		gen dots_binid = _n
		merge 1:1 dots_binid using "temp.dta"
		drop _merge
		erase "temp.dta"

		* Semiparametric Robinson estimate
		semipar mmargin_nbnn_so $controls_educ, ///
			nonpar(trep_date_ulttrans_norm) nograph generate(yhat) ///
			partial(yhat_parout1)

		* Optimal bins for the partialled-out DV 	
		binsreg yhat_parout1 trep_date_ulttrans_norm, ///
			savedata(temp) replace

		* Rename to make room for the new bin locations and binned means
		drop dots_isknot
		rename dots_x dots_x_nocontrols
		rename dots_fit dots_fit_nocontrols

		* Merge new bin locations and binned means
		merge 1:1 dots_binid using "temp.dta"
		drop _merge
		erase "temp.dta"

		* lp robust: partialed out dv
		lprobust yhat_parout1 trep_date_ulttrans_norm, ///
			genvars eval(dots_x)

		renvars lprobust_eval - lprobust_CI_r_rb, postfix(_controls)	
			
		* lp robust: no controls	
		lprobust mmargin_nbnn_so trep_date_ulttrans_norm, ///
			genvars eval(dots_x_nocontrols)

		* X-axis locations for labels
		quietly sum trep_date_ulttrans
		local global_mean = r(mean)
		local global_sd = r(sd)
		local x1 = (tc(20oct2019 17:30:00) - `global_mean') / `global_sd'
		local x2 = (tc(20oct2019 18:30:00) - `global_mean') / `global_sd'
		local x3 = (tc(20oct2019 19:30:00) - `global_mean') / `global_sd'
		local x4 = (tc(20oct2019 20:30:00) - `global_mean') / `global_sd'

		* Create figure
		twoway (scatter dots_fit_nocontrols dots_x_nocontrols, ///
				msymbol(Th) msize(vsmall) mcolor(midblue)) ///
			(line lprobust_gx_bc lprobust_eval, ///
				lcolor(midblue) lwidth(medthick)) ///
			(scatter dots_fit dots_x, ///
				msize(vsmall) mcolor(gs6)) ///
			(line lprobust_gx_bc_controls lprobust_eval_controls, ///
				lcolor(gs6) lwidth(medthick)) ///
			(line lprobust_CI_l_rb lprobust_eval, ///
				lpattern(dash) lcolor(midblue) lwidth(thin)) ///
			(line lprobust_CI_r_rb lprobust_eval, ///
				lpattern(dash) lcolor(midblue) lwidth(thin)) ///
			(line lprobust_CI_l_rb_controls lprobust_eval_controls, ///
				lpattern(dash) lcolor(gs6) lwidth(thin)) ///
			(line lprobust_CI_r_rb_controls lprobust_eval_controls, ///
				lpattern(dash) lcolor(gs6) lwidth(thin)), ///
			legend(order(2 4) ring(0) pos(11) cols(1) lab(2 "No controls") ///
				lab(4 "Education, region, and rurality controls")) /// 
			graphregion(color(white)) ///
			xsize(7) ysize(4.7) ///
			ylabel(-0.1(0.1)0.4, angle(0) nogrid)  ///
			xlabel(`x1' `" "10/20" "5:30 p.m." "' ///
			`x2' `" "10/20" "6:30 p.m." "' ///
			`x3' `" "10/20" "7:30 p.m." "' ///
			`x4' `" "10/20" "8:30 p.m." "', labsize(small)) ///
			ytitle("Average Incumbent Margin") ///
			xtitle("Transmission Time") ///
			title("")
		graph export "outputs/generated/Figure_1a.pdf", replace

	* Figure 1b
	*----------
		use "data/outputs/computo_ulttranstrepdate_missingsatmunmedian.dta", clear	

		* Nulos share
		gen nulosshare = Nulos / total_so

		* Margin
		gen margin_nbnn = mshare_nbnn - cshare_nbnn

		* Recinto ID
		egen reci_id = group(Pais Dep Prov Muni Loc Reci)

		* De-mean
		foreach var of varlist mshare_nbnn cshare_nbnn margin_nbnn nulosshare {
			egen reci_mean = mean(`var'), by(reci_id)
			gen `var'_dm = `var' - reci_mean 
			drop reci_mean
			qui su `var'
			gen `var'_dm_plus = `var'_dm + `r(mean)' // Add overall mean for scale
		}

		* Trimming
		cumul trep_date_ulttrans, gen(cum_date)
		drop if cum_date <= 0.02 | cum_date >= 0.98
		drop cum_date
		
		* Without this, Stata can't calculate the right ROT bandwidth
		su trep_date_ulttrans, d
		local global_mean = r(mean)
		local global_sd = r(sd)
		gen trep_date_ulttrans_norm = (trep_date_ulttrans - `global_mean') / `global_sd'

		* Binned data using binsreg and all booths
		tempfile temp
		binsreg margin_nbnn_dm trep_date_ulttrans_norm, savedata(`temp') 

		gen dots_binid = _n
		merge 1:1 dots_binid using `temp'
		drop _merge

		* nonparametric fit
		lprobust margin_nbnn_dm trep_date_ulttrans_norm, genvars

		* X-axis locations for labels
		sum trep_date_ulttrans
		local global_mean = r(mean)
		local global_sd = r(sd)
		local x1 = (tc(20oct2019 17:30:00) - `global_mean') / `global_sd'
		local x2 = (tc(20oct2019 18:30:00) - `global_mean') / `global_sd'
		local x3 = (tc(20oct2019 19:30:00) - `global_mean') / `global_sd'
		local x4 = (tc(20oct2019 20:30:00) - `global_mean') / `global_sd'
	
		* Create figures
		twoway (scatter dots_fit dots_x, msize(small) mcolor(midblue)) ///
			(line lprobust_gx_bc lprobust_eval, ///
				lcolor(gs6) lwidth(medthick)) ///
			(line lprobust_CI_l_rb lprobust_eval, ///
				lcolor(gs6) lwidth(thin) lpattern(dash)) ///
			(line lprobust_CI_r_rb lprobust_eval, ///
				lcolor(gs6) lwidth(thin) lpattern(dash)), ///
			xlabel(`x1' `" "10/20" "5:30 p.m." "' ///
			`x2' `" "10/20" "6:30 p.m." "' ///
			`x3' `" "10/20" "7:30 p.m." "' ///
			`x4' `" "10/20" "8:30 p.m." "', labsize(small)) ///
			graphregion(color(white)) ///
			ylab(, glcolor(none) angle(0)) ///
			ytitle("De-meaned Average Incumbent Margin") ///
			xtitle("Transmission Time") ///
			xsize(7) ysize(4.7) ///
			legend(off)
		graph export "outputs/generated/Figure_1b.pdf", replace

*------------------------------------------------------------------------------*
timer off 2
timer list 2

log close
clear all

*------------------------------------------------------------------------------*
/* Check that user-written command SEMIPAR works correctly

use http://fmwww.bc.edu/ec-p/data/wooldridge/hprice3, clear

* Command as we are using it in our code
semipar lprice ldist larea lland rooms baths age, ///
		nonpar(linst) nograph partial(yhat_partout_command)

* What the command should be doing, approximately
foreach var of varlist lprice ldist larea lland rooms baths age {
	lpoly `var' linst, at(linst) gen(`var'_hat) degree(1) nograph
	gen `var'_tilde = `var' - `var'_hat
	}
reg lprice_tilde ldist_tilde larea_tilde lland_tilde rooms_tilde baths_tilde age_tilde
gen yhat_partout_byhand = lprice - (_cons + ldist*_b[ldist_tilde] ///
								  + larea*_b[larea_tilde] ///
								  + lland*_b[lland_tilde] ///
								  + rooms*_b[rooms_tilde] ///
								  + baths*_b[baths_tilde] ///
								  + age*_b[age_tilde])

* Compare command output with our version
twoway (lpoly yhat_partout_command linst) (lpoly yhat_partout_byhand linst, yaxis(2)) 
						  
								  






